Dynamics of liquid silica as explained by properties of the potential energy landscape 
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The dynamics of silica displays an Arrhenius temperature dependence, classifying silica as a strong 
glass-former. Using recently developed concepts to analyse the potential energy landscape one can 
get a far-reaching understanding of the long-range transport of silica. It can be expressed in terms 
of properties of the thermodynamics as well as local relaxation processes, thereby extending the 
phenomenological standard picture of a strong glass-former. The local relaxation processes are 
characterized by complex correlated sequences of bond breaking and reformation processes. 

PACS numbers: 64.70.Pf, 65.40.Gr, 66.20.+d 



Silica is a prototypical and technologically relevant 
glass-former, displaying a variety of remarkable physi- 
cal properties like thermodynamic anomalies [lj, |2|, |3|] . In 
contrast to most other glass-formers the temperature de- 
pendence of its transport properties like the oxygen self- 
diffusion constant D(T) display a simple Arrhenius be- 
havior with an activation energy Vdif fusion = 4.7 eV Q 
and is thus a strong glass- former |£||fg. This suggests that 
the transport can be described as a successive breaking 
and reformation processes of Si-0 bonds with an activa- 
tion energy close to Vdif fusion 00- 

To scrutinise this simple picture and thus to obtain a 
microscopic picture of the dynamics of silica we employ 
the framework of the potential energy landscape (PEL), 
defined in the high-dimensional configuration space [8(. 
At low temperatures the properties of silica and other 
glass forming systems are mainly characterized by the 
properties of the local potential energy minima of the 
PEL (denoted inherent structures, IS) [1,0. The ther- 
modynamics of the system is mainly governed by the 
energy distribution of the number of IS. Introducing 
the configurational entropy S C (T) as a measure for the 
number of accessible IS at a given temperature, there 
is an empirical connection of S C (T) to the dynamics 
D(T) oc exp(-A/TS c (T)) (Adam-Gibbs relation [Hj) 
with some fitting parameter A 0, ITTl IT^ . Its theoret- 
ical foundation, however, is under debate[l3L Il4| and 
no direct interpretation of A is available, yet. In any 
event, one would expect that also the topology of the 
PEL should be of utmost importance for understanding 
the dynamics [l5j . 

Following the ideas of Stillinger and Weber one can 
map a trajectory, obtained via a molecular dynamics sim- 
ulation, to a sequence of IS via frequent minimization 
of the potential energy. One example is shown in Fig.l. 
This highlights how the system is exploring the PEL. One 
can group together IS to metabasins (MB) such that the 
dynamics between MB does not contain any correlated 
backward- forward processes [1,E1>E1- Each MB is char- 
acterized by a waiting time r and an energy e (defined as 
the lowest IS energy in this MB). The technical details of 
this approach, which have been developed in the context 
of studying a simple glass-forming model, i.e. the binary 
mixture Lennard- Jones system (BMLJ), can be found in 
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Figure 1: Time dependence of the IS energy e during a molec- 
ular dynamics simulation at T = 3000 K. The fountain-like 
objects correspond to time periods during which the system 
is moving very fast in configuration space. e c is an estimate 
for the low-energy cutoff of the PEL. 



In this work we apply these techniques to silica. Com- 
bining information about the energy distribution of IS 
and the local relaxation processes, reflecting the local 
topology of the PEL, we obtain a far-reaching under- 
standing of its dynamics. From this we can identify the 
reasons why silica is a strong system and obtain a quan- 
titative understanding of observables like the resulting 
activation energy Vdif fusion- Further information is ob- 
tained from an appropriate comparison with the real- 
space behavior of silica. The underlying picture, emerg- 
ing from these results, extends substantially the rational- 
ization of the strong behavior of silica, sketched above. 

Information about the BKS potential [l8|. used to 
model silica, as well as further simulation details can be 
found in Ref. 0|. For an optimum analysis in terms of 
the PEL the system size should be as small as possible 
without showing significant finite size effects ^(|. It has 
been shown that already for system sizes N ss 10 2 finite 
size effects concerning the configurational entropy [l9j . 
the properties of tunneling systems I20| , the temperature 
dependence of the oxygen diffusion [ijj as well as the na- 
ture of the relaxation processes in BKS silica (checked, 
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Figure 2: Temperature dependence of the oxygen diffusion 
constant D(T) (macroscopic transport) and the inverse av- 
erage waiting time (r(T)) -1 (microscopic relaxation). The 
low-temperature activation energy is Vdif fusion = 4.84 eV. 
Around 3500 K the so-called fragile-to-strong crossover is ob- 
served [j sum!. 



e.g., via the degree of non-exponentiality in the incoher- 
ent scattering function) |19ll2l] |are small in the accessible 
range of temperatures. Here we choose N = 99. Proper- 
ties of larger systems can be then predicted from statis- 
tical arguments 22] . Recent studies have shown that the 
distribution of configurational energies has a low-energy 
cutoff around some energy e c with a finite configura- 
tional entropy Q . It results from the network constraints 
in defect-free configurations It will turn out to be 

one key feature for the understanding of the dynamics. 

In analogy to previous results for the BMLJ system 
[23l ] the oxygen diffusion constant D(T) is proportional 
to the inverse of the average MB waiting time (r(T)); 
see Fig. 2. Thus, a local quantity like (r(T)} fully de- 
termines the temperature dependence of diffusion, i.e. 
d\nD{T) / d(3 . The low-temperature activation energy 
Vdif fusion — 4.84 eV is very close to the experimentally 
observed value of 4.7 eV 0. Around 3500 K one observes 
the crossover from the high-temperature non-Arrhenius 
to the low-temperature Arrhenius-regime Q, |jj [T2I |24| . 

One may suspect that the MB waiting times are 
strongly energy dependent because high-energy config- 
urations have a larger number of defects which can 
more easily relax |20|. For a quantification we deter- 
mine the average MB waiting time in dependence of tem- 
perature and energy, denoted (r(e,T)) from analysis of 
several independent equilibrium runs; see Fig. 3. Inter- 
estingly, for all e one finds a simple Arrhenius-behavior, 
characterized by an effective activation energy Vmb{z) 
and a prefactor Tn(e). Thus, the Arrhenius-behavior 
is not only present for simple atomic systems like the 
BMLJ system [ll| but also for silica. One observes a 
crossover from the high-energy regime with Vmb{&) ~ 
Vq and To(e) = T m i cro to a low-energy regime with a 
strong energy-dependence of both functions. The result- 
ing broad distribution of waiting times implies that en- 
ergy is likely the most dominating factor for understand- 
ing the occurrence of dynamic heterogeneities in silica 
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Figure 3: Arrhenius plot of the average MB waiting time 
(r(e, T)) in dependence of energy (the chosen energy bins can 
be read off from the inset). (r(e,T)} can be characterized by 
an effective activation energy Vmb(&) and a prefactor ro(e), 
shown in the inset. There exists a dynamic crossover energy 
e cross such that for e > e cross , r (e) = r micro ?s 20 fs and 
VifB(e) = Vb ~ 0.8 eV are roughly constant. The lines are 
guides to the eyes. 



and other glass- forming systems [25j . 

For relating the thermodynamics and the dynamics we 
introduce p(e, T) as the Boltzmann probability to be in a 
MB of energy e. It turns out that in the relevant energy 
and temperature range p(e,T) is virtually identical for 
MB and IS 0] and thus contains all the information 
about the configurational entropy. Using the quantities, 
introduced so far, we can write down the formal relation 
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Its physical relevance is far-reaching because it relates the 
thermodynamics (viap(e, T)) and the local dynamics (via 
(r(e,T))) to the long range diffusion. It follows that for 
very low temperatures for which p(e,T) has only contri- 
butions for e w e c the dominant contributions to the av- 
erage waiting time originate from configurations with en- 
ergies close to e c via (r(e « e c , T)). Then the local Arrhe- 
nius behavior (r(e « e c , T)) cx exp(/3VMs(e w e c )) trans- 
lates into a macroscopic Arrhenius behavior. Indeed, 
the macroscopic activation energy Vdif fusion is close to 
V M B(e w e c ); see Fig. 4. 

What determines the value of the crossover tempera- 
ture of 3500 K? At this temperature p(e,T) is peaked 
around e c + 4 eV and the low-energy wing of p(e, T) 
just starts to be influenced by the low-energy cutoff |19j . 
Thus, on first view the above arguments to rational- 
ize Arrhenius behavior should only apply for lower tem- 
peratures. However, due to the additional weighting of 
exp(/3VAfs(e)) by I/r (e) in Eq.l, which for e « e c is 
more than four orders of magnitude larger than I / T m i cro 
(ro(e we c )w5 - 10~ 3 fs), the influence of the low-energy 
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Figure 4: Detailed analysis of the escape properties from 
one defect-free low-energy structure ISo- a) Temperature de- 
pendence of the average waiting time in ISo- b) One specific 
escape path starting from ISo- This representation reflects 
the sequence of inherent structure energies as well as saddle 
energies. The results for pback are given above the arrows. 
The pairs of numbers at the bottom reflect the numbers of 
silicon and oxygen defects, respectively. The resulting effec- 
tive barrier height E e3cap e,a is indicated, c) Distribution of 
Eescape 7 a and its average value (Eescape) — 4.6 eV Vescape) 
for 30 different escape paths from ISo. 

states in the integral is significantly enhanced, thus giving 
rise to the actually observed crossover. Is there a direct 
relation to the mode-coupling temperature T c w 3300A" 
p4|? Qualitatively, T c is related to the beginning domi- 
nance of activated processes rather than to the presence 
of a low-energy cutoff of the PEL and thus seems to have 
a different physical origin. 

In the next step we elucidate the microscopic origin of 
the escape properties from configurations close to e c . We 
present detailed results for one specific configuration (de- 
noted ISo). In a first step we determine its average wait- 
ing time via repeated runs starting at ISo with different 
initial velocities and different temperatures. These runs 
typically involve many unsuccessful escape attempts. We 
define the waiting time as the time when ISo is left for the 
last time. We find a simple Arrhenius behavior with an 
effective activation energy of V esca pe — 4.77 eV (Fig. 4a). 

This activation energy incorporates all the complexity 
of the local PEL around this configuration. In partic- 
ular, the configuration may relax along different paths 
a, which all contribute to the total relaxation rate, i.e. 
r = Y^, -T" • 11 Pa is the probability to escape via path 
a and if path a is characterized by an effective barrier 
height E esca p e>a (see below for its definition) one can 
write (E escape ) = ^2p a E escapeia |l6], which is the aver- 
age value of E escape ^ a for independent runs starting from 
ISo- A typical escape path is shown in Fig. 4b. One can 
see the sequence of IS after ISo is visited for the last time. 
For every ISj a value, denoted Pback, has been obtained 
from counting for a set of independent simulations with 
starting structure IS,, whether or not the system returns 
to ISq. Qualitatively, Pback < 0.5 indicates that the sys- 



tem will (on average) escape from the catchment region 
of ISo- This limit typically involves the entropic effect 
that due to the large number of transition options in 
the high-dimensional PEL there is no need to follow the 
path back to ISo- As shown in Ref. [16j for a given escape 
path a good estimate of E escape .a is first to locate the 
first IS with pback < 0.5 and then to identify E escape ,a as 
the highest energy along the reaction coordinate up to 
this IS. Application to the escape path in Fig. 4b yields 
E e scape,a = 4.3 eV. Four transitions are required until for 
the first time pback < 0.5. 

From the distribution of E escape , a for 30 different es- 
cape paths from ISo one obtains (E escape ) = 4.6 eV which 
is close to the value of V escape = 4.77 eV ; see Fig. 4c. 
Thus, it is indeed possible to quantitatively relate the 
effective activation energy to the local properties of the 
PEL. Repeating this analysis for two different low-energy 
IS we get a similar agreement. More generally this im- 
plies that Vmb(s ~ e c ), on the one hand, can be quan- 
titatively related to specific local barriers and, on the 
other hand, to the activation energy of macroscopic dif- 
fusion Vdif fusion- This establishes a direct link between 
the microscopic and macroscopic behavior of silica. 

Having identified a complete escape process (relevant 
for Vdif fusion) via the pfc ac fe-criterion we are in a posi- 
tion also to analyze its real-space characteristics. First 
we note that in the vast majority of cases the sequence 
of IS-transitions during the escape is correlated. This is 
reflected by the fact that at least one atom, involved in 
a bond-breaking or reformation process during a specific 
IS-transition (i > 1) ISj — >IS;+i was involved during a 
previous IS-transition. When comparing IS4 with ISo hi 
total 4 Si-0 bonds have been broken and 5 silicon atoms 
have changed their oxygen neighbors. These values are 
close to the average behavior after analyzing the escape 
from all three three initial IS (4.6 and 4.4, resp.). This 
implies significant correlated bond-breaking and reforma- 
tion processes. In particular, Vdif fusion cannot be related 
to the breaking of a single Si-0 bond. Other researchers 
have rationalized the value of Vdif fusion by the sum of 
half of the mean formation e nerg y of an oxygen Frenkel 
pair and a migration barrier 21] . On a qualitative level 
something similar is observed here because in the first 
step a defect is created and afterwards the defect is trans- 
ferred until pback < 0.5. A closer look, however, reveals 
that the behavior in BKS silica is more complex as re- 
flected, e.g., by the fact that in the example of Fig. 4 in 
between also configurations with no defects occur (IS3). 

The additional contribution of the saddle between IS3 
and IS4 of approx. 0.8 eV to the final value of {E escape ) is 
small. This also holds in general (approx. 1.0 eV) which, 
interestingly, is close to Vq . Thus one may conclude that 
there are two distinct contributions to the activation en- 
ergy Vmb{&)'- (1) VMs(e) — Vb as the contribution re- 
flecting the topology of the PEL, related to differences 
between IS energies, and (2) Vq as the additional con- 
tribution of the final saddle. Around e ~ e c the first 
contribution is dominating. 
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Why is silica a strong glass-former? The standard sce- 
nario, sketched in the introductory paragraph, would im- 
Ply (r(e,T)) 

f fusion 

) lalg. This re- 
flects the presence of one typical relaxation process which 
holds throughout the entire PEL [8J. In contrast, our 
simulations have revealed a strong energy dependence 
of Vmb(g) and ro(e) together with complex successive 
bond-breaking and reformation processes. Rather we can 
identify two underlying reasons for the classification of 
silica as a strong glass-former: (1) the presence of the 
cutoff in the PEL of silica as a consequence of its net- 
work structure, (2) the Arrhenius temperature depen- 
dence of (r(e k, e c ,T)) together with a large attempt 
rate l/ro(e c ). How to rationalize property (2)? As can 
be seen from Fig. 4c the distribution of effective barrier 
heights E escape a is very narrow. For the example of 
ISo this implies an Arrhenius temperature dependence 
in agreement with the observation; see Fig. 4a. This be- 
havior was also observed for the escape from the other 
low-energy IS, analyzed along the same lines. Thus, it 
seems to be a general feature that starting from a low- 
energy (and typically defect-free) configuration around 
e c the system first has to acquire an energy of approx. 
e ~ e c + Vmb{g c ) until the escape is complete (ending in 
an IS with e w e c + Vmb{e c ) — Vq)- More pictorially one 
may state that low-energy IS (e « e c ) form the bottom 
of crater-like objects in the PEL and the system has no 
escape option apart from climbing up the whole crater 
until e rs e c + Vmb{zc)- Thus, beyond the presence of 
the low-energy bound of the PEL the strong behavior of 
silica is also related to this crater-like structure of the 
PEL. 

The previous thermodynamic analysis has revealed 
that for this system the number of IS with e c + Vmb (ec) — 
Vo w e c + 4 eV is approx. 10 5 times larger than the cor- 
responding number at e rj e c 0] . This observation sug- 
gests that the number of possible transitions from e « e c 
to configurations with energies around e c + 4 eV is also 
exponentially large, thereby rationalizing the dramatic 
increase of l/ro(e c ) as compared to l/T m i cro . It has been 
already explicitly shown in previous model calculations 
that the prefactor for relaxation processes in model land- 
scapes scales with the number of accessible states for 
which for the first time Pback is smaller than approx. 50% 



[26j| . Then l/ro(e c ) can be much larger than microscopic 
jump rates. Interestingly, in the relevant energy regime 
of the BMLJ system the number of IS increases much 
weaker with increasing energy and correspondingly there 
is hardly any energy-dependence of To(e) [lfj. This fur- 
ther supports our hypothesis about the origin of the large 
energy-dependence of To(e). 

It has been speculated that different network-forming 
systems (e.g. silica and water) have similar properties[27l 
Ea. l29j . Indeed, indications have been found recently 
that also the amorphous states of water possesses a low- 
energy cutoff, which will influence the thermodynamics 
and dynamics at low temperatures similarly to silica j^. 
Therefore exploration of the properties of silica may be 
of major importance also for an improved understanding 
of water and other network formers. The possible uni- 
versality of this class of systems has recently lead to the 
formulation of an abstract model which is aimed to re- 
flect the basic physics of all amorphous network formers 

m 

In summary, using the PEL-framcwork we have iden- 
tified relevant underlying mechanisms for the dynamics 
of silica: (1) The elementary relaxation processes are 
not simple bond-breaking processes with the final acti- 
vation energy but rather form correlated sequences of 
many bond-breaking processes with a strong dependence 
on the initial energy. (2) The simultaneous presence of 
the low-energy cutoff of the PEL as well as the narrow 
distribution of escape barriers from MB give rise to the 
resulting strong behavior. (3) Dramatic entropic effects 
are present for the escape from low-energy states, show- 
ing that the dynamics is much more complex than re- 
flected by the resulting simple Arrhenius behavior. (4) 
The occurrence of the crossover temperature around 3500 
K can be quantitatively understood by the presence of the 
low-energy cutoff and these entropic effects. (5) There ex- 
ists a crossover energy which separates the high-energy 
liquid-like behavior and the low-energy activated behav- 
ior. 
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